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ABSTRACT 

Dispersion in the interstellar medium is a well known phenomenon that follows a simple relationship, which has been used to predict 
the time delay of dispersed radio pulses since the late 1960s. We performed wide-band simultaneous observations of four pulsars 
with LOFAR (at 40-190 MHz), the 76-m Lovell Telescope (at 1400 MHz) and the Effelsberg 100-m Telescope (at 8000 MHz) to 
test the accuracy of the dispersion law over a broad frequency range. In this paper we present the results of these observations which 
show that the dispersion law is accurate to better than 1 part in 10 5 across our observing band. We use this fact to constrain some of 
the properties of the ISM along the line-of-sight and use the lack of any aberration or retardation effects to determine upper limits 
on emission heights in the pulsar magnetosphere. We also discuss the effect of pulse profile evolution on our observations, and the 
implications that it could have for precision pulsar timing projects such as the detection of gravitational waves with pulsar timing 
arrays. 
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1. Introduction 

As radio emission passes through the interstellar medium (ISM) 
it interacts with electrons, which cause it to be dispersed. This 
is observable through pulsed emission, where a low-frequency 
pulse is delayed with respect to the same pulse at higher fre- 
quencies. For a pulse at frequency v (frequency will be given 
in MHz here and throughout the paper) travelling a distance, D, 
through an unmagnetised ionised gas, the dispersive delay with 
respect to infinite frequency, A?dm, is given by: 



( f° dl \ D 

Af D M = — , 

Wo V g ) c 



where v g is the group velocity: 



v p is the plasma frequency: 

(in cgs units) 
8.98 x 1(T 3 V«7MHz, 




(1) 

(2) 

(3) 
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c is the speed of light, m e and e are the charge and mass of 
an electron respectively, and n e is the electron density in cm" 3 . 
Normally, this time delay can be approximated using the first 
term in the Taylor expansion of 1 / v g , giving the cold dispersion 
law3 
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DM 



2.41 x 10- 



(5) 



where DM is the dispersion measure, the integrated column den- 
sity of free electrons in pc cm" 3 . This relation was used as early 
as 1968 to accurately predict the dispersive ti me delay to within 
1 part in 3000 between 40 MHz and 430 MHz ( Tanenbau m et alj 
[T968h . However, as ever higher timing precision is required 
for projects like u sing pulsars for gravitational wave detection 
dJenet et al.l2005l) it is important that any second-order effects of 
the ISM, such as refractive delays, DM variation s and delays as- 
sociated with pulse broad ening from scattering (Foster & Cordesl 



1990; 



Cordes & Shannon 2010), are also studied and understood 



fully (Yo u etal.11200 8; Hemb erger & Stineb ring 2008). Most of 



these proposed effects have strong frequency dependencies, with 
scaling indices between v~ 3 and v~ 4 , and are therefore most 
prominent at low frequencies. 

As well as second-order ISM effects, there are other 
proposed frequency-dependent effects, such as propaga - 
tion from within the pulsar magnetospher e dMichell fl 99 1 ), 
super-dispersion dshitov & Malofeevl fl985l iKuz'minl 11986: 
Shitov et al.ll98^lKuz^rm : ri*eTan2*0*0*8T) and aberration and retar- 
dation dCordesl!978l) . The ionosphere also contributes slightly to 
the DM but typically the total number of electrons alo ng a path 
throu gh the ionosphere is less than 100 TEC units (e.g. lLiu et al.l 

"In 



120091) . which corresponds to a DM of just 3 .24 x 1 0~ 5 pc cm 3 . 
practice, this contribution is small and follows a v~ 2 law, making 
it indistinguishable from the dispersive delay of the ISM. 



1 Note that often, the more precise value of 2.410332 x 10" 4 is used 
as the constant in Equation[5] however in this paper we follow the con- 
vention of using 2.41 x 10~ 4 . This is only important when considering 
absolute values of DM. It should also be noted that this approximation 
only holds when the observing frequency is much greater than both the 
local plasma frequency and electron gyrofrequency. 



As the sizes of these effects have not yet been determined, 
their frequency dependence could introduce systematic errors 
when data from different frequencies are combined. The size of 
these effects could also potentially vary with time. For exam- 
ple, if the distribution of free electrons along the line-of-sight 
changes, the magnitude of some of these effects will vary, intro- 
ducing further errors into pulsar timing data. 

Another potential source of error arises from pulse profile 
evolution. Most (if not all) pulsars show variation in the shap e 
of their pulse profile as a function of frequency (ICraftl [T970b . 
This change can be anything from a slight broadening of the 
pulse, or difference in the relative positions of the components, 
to components appearing or disappearing completely. 

Even though wide-band simultaneous observations 
have bee n possible to perform for a long time (see for 
example iPhilUps & Wolszczanl Il992t IKuz'min et at] 1 19981 : 
lHankins & Rankinll2010l) the relatively narrow frequency bands 
compared to the separation interval between frequencies have 
meant identifying exactly how the components evolve has been 
difficult, particularly at the lowest frequencies. One of the 
reasons for this is that profile evolution, which only manifests 
itself clearly over very wide frequency ranges, is difficult to 
disentangle from dispersion (dispersion is normally removed by 
fitting the data with a v~ 2 power l a w wh ich mimics Equation 
|5). For example, in IKuz'min et alj (1998) the authors suggest 
that PSR B0809+74 appears to have an extra, non-dispersive 
delay of ~ 30 ms at 10 GHz, though this could alternatively be 
explained by pulse profile evolution. 

To calculate pulse times of arrival (TOAs), the data are cross- 
correlated with a 'template', a smoothed model of the pulse pro- 
file. The peak of this cross-correlation spectrum is used to de- 
termine when the pulse arrives relative to the template, and the 
phase of the Fourier transform of the cross-correlation is then 
used to measure the phase offset. If the pulse profile evolves sig- 
nificantly with frequency the shape of the data and the template 
can be slightly different. These subtle changes in shape can in- 
troduce frequency-d ependent errors into TOAs if they are not ac- 
counted for properly. lAhuia et al.l(l2007l) simulated the effects of 
pulse profile evolution on measurements of DM. They found that 
using a template which is slightly different from the shape of the 
data can cause a gradient in the phase of the cross-correlation, 
which causes an offset in the TOA. 

Some of the simulated errors, in the work of Ahuja et al, due 
to this effect in normal (slow) pulsars were as large as millisec- 
onds, although the magnitude of the effect on real data has not 
been studied in detail until the work presented here. The size of 
this effect should scale linearly with pulse period, and so it is ex- 
pected to be much smaller for millisecond pulsars, which is why 
sub-micro second timing prec ision is already common (for ex- 
ample, see lVerbiest et al.l2009l) . But, with the introduction of the 
first ultra-broadband receivers to be u sed for pulsar timing (e.g. 
iDuPlain et alj r2008; Jones et al. 2010), which will be capable of 
observing with large fractional bandwidths even at frequencies 
of a few gigahertz, pulse shape variation across the band will 
become far more pronounced, and this effect will become more 
apparent. 

We took observations of four bright pulsars (pulsars 
B0329+54, B0809+74, Bl 133+16 and B1919+21) simultane- 
ously at multiple frequencies with the LOw Frequency ARray 
(LOFAR, see van Haarlem et al. in prep.), the 76-m Lovell 
Telescope and the Effelsberg 100-m Telescope, spanning a fre- 
quency range between 40 and 8350 MHz to try to constrain the 
properties of the ISM and test the cold plasma dispersion law. 
Details of these observations are given in Section [2] Our obser- 
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vations cover a very large frequency range simultaneously and 
each recorded band has a significant fractional bandwidth, also 
making them useful for studying pulse profile evolution. LOFAR 
provides the largest fractional bandwidth which has ever been 
possible at the lowest radio frequencies observable from Earth. 

The observations from both LOFAR observing bands were 
also taken truly simultaneously, both starting and ending at ex- 
actly the same tim^E and passing through the same hardware so 
no timing model was required to align the profiles precisely (as 
opposed to using overlapping observations with a timing model). 
This part of our frequency range provides the best test of steeply 
scaling frequency-dependent delays and this allows us to place 
strong constraints on the magnitude of the second-order ISM ef- 
fects. Any deviation from the v~ 2 law larger than a few millisec- 
onds would be clearly visible over this large bandwidth and at 
these low frequencies. 

Similar observations have been carried out previously by 
iTanenbaum et alJ dl968l) . who were able to set a surprisingly 
good limit that the error in the dispersion equation, At < 3 x 
1(T 4 s between 40 MHz and 430 MHz just six months af- 
ter the discovery o f pulsars in 1967. Since t hen similar work 
has been done by Shitov & Mal ofeevl (1 19851) . who found evi- 
dence for exc ess dispersion at low frequencies. This was also 
supported by lHankins et al.l (Il99ll) . who noticed a difference 
in the DM values determined from aligning single pulses and 
those determined from aligning the average pulse profile. These 
delays at low frequencies have been termed 'super-dispersion' 
and attri buted to m agnetic sweepback i n the pulsar magneto- 
sphere (IShitovll 983). Subsequent work by Phillips & Wols zczanl 
d!992h and lAhuia et al.l (120051) found contradictory results, with 
no evidence for extra dispersion at low frequencies, but some 
pulsars (B0525+21, B 1642-03 and B 1237+25) showed an in- 
creased dispersion measure at high frequencies. These high fre- 
quency delays were explained by propagation delays in the pul- 
sar magnetosphere. 

In this paper, after detailing the observations (Section|2}, our 
analysis (Section [3), and the simulations used to determine up- 
per limits on the magnitude of the delay present in our data 
(Section |4|; we determine the effect of the ISM on pulsar tim- 
ing and extract information about the composition of the ISM 
(Section [5]) and the pulsar magnetosphere (Section |SJ. We also 
use these data to determine how much of an impact profile evo- 
lution has on pulsar timing and use the pulse profiles in each of 
the frequency bands to construct a model of how the pulse pro- 
file evolves with frequency (Section |7). Finally, we summarise 
our main findings in Section[8] 

2. Observations 

On 1 1 December 2009 we observed four pulsars using LOFAR, 
the 76-m Love 11 Telescope and the Effelsberg 100-m Telescope 
simultaneously. LOFAR is an international interferometric tele- 
scope, comprised of many thousands of dipole antennas grouped 
into 'stations' and operating in the lowest four octaves of the 'ra- 
dio window' visible from the Earth's surface. The LOFAR sta- 
tions are arranged in a sparse array, spread across Europe, with 
a dense core region located in the Netherlands. LOFAR operates 
in two frequency bands which straddle the FM radio band. The 
Low Band Antennas (LBAs) can record 48 MHz of bandwidth 
between 10 and 90 MHz and the High Band Antennas (HBAs) 
can record 48 MHz of bandwidth between 110 and 240 MHz. 

2 There are small clock offsets between the LOFAR stations at this 
stage in construction, but these are less than 20 ns. 



Table 1. Data characteristics. Shown are the centre frequency 
(v), bandwidth (B), number of channels {N c h- dn ), and the sampling 
time (r samp ) for the data of each of the four telescopes used in 
the simultaneous observations. 





v (MHz) 


B (MHz) 


chan 


Fsanip (ms) 


DE601 


60.15625 


36.328125 


2976 


1.31072 


CS302 


163.28125 


48.4375 


3968 


1.31072 


Lovell 


1524.0 


512 


512 


1.0 


Effelsberg 


8350.0 


1000 


1° 


~ 1.0* 



Notes. (<!) Data was written as a single dedispersed timeseries. This 
varied depending on the source. The pulse period was divided into 1024 
bins. 

Table 2. Pulsar characteristics. Given are the integ ration time 
(Tint) , catalogue value for dispersion measure (DM, Hobbs ~et~al1 
2004), the pulse period (P, determined from regular tim- 
ing observations from the Lovell telescope at Jodrell Bank 
Observatory), the LOFAR observation ID and the number of bins 
across the pulse profile (AW for each of the pulsars observed.) 



Tint DM P LOFAR N° m 

(s) (pc cm' 3 ) (s) Obs ID 

B0329+54 7200 26.833 0.714536 L2009_16116 256 

B0809+74 5400 6.116 1.292209 L2009_16102 512 

B1133+16 10800 4.864 1.187799 L2009_16100 512 

B1919+21 5400 12.455 1.337360 L2009_16104 512 



Notes. (a) Although the pulsars in our sample were observed with higher 
time resolution, the data were downsampled so that each observation 
had the same number of bins in the pulse profile, so that the data from 
different frequencies could be compared directly. 



For a full description of how LO FAR is used for pulsar obser- 
vations see Stappers et al] d201 ll) and for a general LOFAR de- 
scription see van Haarlem et al. (in prep.). 

At the time of the observations, LOFAR was still under con- 
struction so only single stations were used to record data and 
the LBAs were only able to record 36 MHz of bandwidth. A 
core station (CS302, which consists of 48 HBA tiles) was used 
to observe in the LOFAR high band between 138.9709 and 
187.4084 MHz and an international station (DE601, based at the 
Effelsberg Radio Observatory and consisting of 96 LBAs) was 
used to observe between 42.0959 and 78.4119 MHz in the low 
band. It is worth noting that the LOFAR stations which were 
used were not yet internally calibrated and hence the sensitivity 
of the observations presented here is significantly lower com- 
pared to what is now possible. 

Data were taken simultaneously with the 76-m Lovell 
Telescope, using the digital filterbank backend in search mode, 
at a central frequency of 1524 MHz with 512 MHz bandwidth 
(see lHobbs et al.ll2004l for details). Simultaneous data were also 
taken with the Effelsberg 100-m Telescope, using the Effelsberg - 
Berkeley Pulsar Processor to coherently dedisperse the data on- 
line, at a c entral frequency of 8350 MHz with 1000 MHz band- 
width (see lBacker et al.lfl997l for details). The observational pa- 
rameters are summarised in Tables Q~]and|2] 

In this paper, we concentrate mainly on the total inten- 
sity profiles of pulsars, and only Stokes I data were recorded 
with each instrument (LOFAR, the Lovell telescope and the 
Effelsberg telescope). However, if there is a significant gain dif- 
ference between the two hands of (either linear or circular) po- 
larisation, or leakage, and the pulsar is strongly polarised then 
there may be distortions in the profile. Instrumental polarisa- 
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Fig. 1. Comparison of timing residuals (TOAs subtracted from a model of the pulsar's expected TOAs) obtained using a single 
template for each frequency band (left), and using a frequency-dependent template (right), plotted against v on a logarithmic scale. 
The residuals from the static templates show significant deviations from white noise. These systematic errors are mostly removed 
with the frequency-dependent templates and the residuals appear straighter, and agree on a single value of DM (apart from those of 
PSR Bl 133+16, see text for details). 'Gaps' in the frequency coverage are from subbands which have been removed because they 
contain strong RFI. The pulsars used are (from top to bottom) B0329+54, B0809+74, Bl 133+16 and B1919+21. 



tion from the Lovell and Effelsbe rg telescopes is known to be 
small (see e.g. iMontene gro-Mon tes et"al] 120081 : iGould & Lvnd 
Il998l) . but at the time of these observations, the LOFAR polari- 
sation data was uncalibrated. We note however that in the case of 
LOFAR the leakage terms are at worst a few percent, and this is 
further reduced due to the fact that there are many thousands of 
elements, which are physically identical, so any leakages would 
average out over the array. Also, as the orientation of the LOFAR 
dipoles is at 45 degrees to the north-south orientation and our 
sources were all observed close to transit, both sets of dipoles 
would have received approximately the same amount of radi- 
ation from the source. Lastly, we note that the polarised frac- 



tion is relatively low, and that the Faraday rotation in the ISM is 
large at these frequencies. In general, with the bandwidths used 
to make pulse profiles (typically ~ 12 MHz, see Section [3]), the 
plane of polarisation is rotated through at least 180°in each pro- 
file, thereby further reducing the effect of polarisation calibration 
terms. We are therefore confident that the results presented here 
are not affected significantly by polarisation calibration. 

3. Analysis 

The data were converted into PSRFI TS format and pro cessed us- 
ing the PSRCHIVE software suite dHotan et alJl2004h . Regular 
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timing observ ations from the Lo vell telescope at Jodrell Bank 
Observatory (Hobbs et aD l2004T) were used to derive accurate 
values for the spin period (P) and spindown (P) for the day of our 
wide -band observ i ng cam paign. Initial DM estimates were taken 
from Hobbs et al.l d2004) (see Table [2]for the values used). The 
data were dedispersed, folded and terrestrial interference signals 
were removed by hand using the RFI-removal software, paz0. 

Templates were made for the observations at each telescope 
by completely collapsing the data in both time and frequency to 
make an average pulse profile. The profiles were then fitted with 
von Mises function^ to create analytic templates using paas 3 . 

A template was created for each observing band and these 
templates were subsequently aligned by eye with pas 3 using ei- 
ther the brightest peak (for PSR B0329+54) or the midpoint 
between the two brightest components (for pulsars B0809+74, 
Bl 133+16 a nd B1919+21) as the fiducial point as described in 
ICraftl JT970b . The templates were cross-correlated with the data 
to get times of arrival (TOAs) using pat 3 and barycentred using 
TEMPofl These TOAs were subtracted from a model of the pul- 
sar's rotation in order to produce 'timing residuals', which are 
plotted in the left hand panel of FigureQ] 

Currently, it is not possible to perform absolute timing with 
LOFAR. The clock corrections between the stations are known 
to better than a few nanoseconds, but the offset between the 
LOFAR clocks and Coordinated Universal Time (UTC) is not 
logged. This introduces time delays which are on the order of 
a few milliseconds between LOFAR data and the data from 
the other telescopes in the observations. The delays presently 
require an arbitrary phase offset ('jump') to be removed be- 
tween data sets from different telescopes. In the case of the 
8.35 GHz Effelsberg observation, where no frequency resolution 
was stored (see Table [TJ, the determination of this phase offset 
uniquely defines the residual with respect to the other observa- 
tions, implying that no further timing information can be derived 
from the observation. We have therefore omitted the Effelsberg 
observations from the timing analysis (though they were still 
used in the analysis of profile evolution, see Section|7]i. 

One can see in Figure Q] (left) that for each of the pulsars 
there is clear structure in the timing residuals as a function of 
frequency. Often, this structure shows different slopes for each 
of the bands, indicating that the deviations from a good fit are 
not simply caused by an incorrect dispersion measure. Initially it 
seemed likely this might be caused by deviations from the simple 
form of the cold plasma dispersion relation. 

Consequently, we tried fitting the data using power laws with 
exponents between -5 and +5, but none of the fits significantly 
improved upon the chi-squared obtained from only fitting for the 
dispersive delay (for example, in the case of PSR B 1 9 1 9 +2 1 , the 
best power-law fit only decreased the reduced chi-squared from 
1 .47 to 1 .44, which, given the number of degrees of freedom in 
the model, still corresponds to a ~ 99% chance that there is un- 
modelled structure in the residuals). In some cases, the residuals 

3 A tool from the PSRCHIVE suite dHotan et alj 12004) . 
http://psrchive.sourceforge.net/ 

4 A von Mises function is given by f(x) = ^fjfi > where /i and 1 Ik 
are analogous to the mean and the variance in a normal distribution. Iq 
is the modified Bessel function. They are used by paas because they are 
needed to deal with pulsars which have broad pulse profiles, although 
for the pulsars in our sample, which have narrow pulse profiles, (the 
pulse width in all cases is < 20% of the pulse period) they are almost 
identical to Gaussians. 

5 tempo2 is a pulsar timing p ackage for b arycen - 
tring and modelling TOAs (iHobbsetalJ [2006). 
http://www.atnf.csiro.au/research/pulsar/tempo2 




Phase Lag O 

Fig. 2. Errors in timing PSR B 19 19+21. The top panel shows the 
HBA template (black line) and the 145 MHz data (grey line). The 
bottom panel shows the cross-spectrum phase of the template 
and the profile and a straight line fit to the data (solid line). The 
subtle difference between the shape of the pulse profile and the 
template causes a gradient in the cross-correlation phase shifting 
the apparent TOA. 



also showed structure which could not be explained by a sin- 
gle function (for example, a positive gradient in the LBA data, 
a negative gradient in the HBA data and a positive gradient in 
the data from Jodrell Bank). No simple power law can account 
for the different slopes in our timing residuals, so we conclude 
that, the structure in the residuals is not caused by the ISM or 
aberration and retardation in the pulsar magnetosphere. 

The real cause of the systematic errors is the p ulse profile 
chang ing as a function of frequency, as noted in lAhuia etaD 
(2007). As the profile changes, the template which is used in the 
cross-correlation becomes more and more inaccurate and there is 
a systematic drift in the residuals towards earlier or later TOAs. 
Figure [2] shows how these errors arise. The top panel shows the 
template used for the HBA observation of PSR B 1 9 1 9+2 1 along 
with the observed pulse profile at 145 MH z. The botto m panel 
shows the cross-spectrum phase (see lAhuia etal]|2007l for de- 
tails) of the two profiles and a straight line fitted to the data. 
There is a gradient in the fitted line which means that the peak 
of the amplitude of the cross-correlation (used to produce the 
TOAs) is shifted slightly, causing the apparent TOA to be de- 
layed or advanced compared to its true value. 

This eff e ct has only been brought to attention recently by 
lAhuia et alJ (120071) . because pulsar timing is normally done at 
high frequencies (where the pulse shape changes are less ob- 
vious), and with modest bandwidths. With these narrow band- 
widths, it is easy to incorrectly interpret this structure as an in- 
correct DM. With the wide fractional bandwidth and high fre- 
quency resolution of our simultaneous data, however, the effect 
cannot be fitted with a simple v~ 2 law, and so the structure is 
more easily identified. 

To remove the effects of profile evolution and find the true 
DM value and residuals, we first modelled how the pulse pro- 
files evolve as a function of frequency. The data were divided 
into narrower-frequency subbands. Each subband was chosen 
such that the profile evolution within it was not significant, but 
there was still enough bandwidth so that the pulsar was detected 
clearly. The size of these bands depended on the signal-to-noise 
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Table 3. Apparent DM in each frequency band of our observations using a static template and a template based on a frequency- 
dependent model. DM avg is the DM using all of the observing bands. 







DM LBA 


DMhba 


DM Jod 


DM 


avg 


B0329+54 


static 


26.774 i 


: 0.001 


26.664 ± 


0.0003 


26.7 


±0.2 


26.7673 4 


: 0.0002 




model 


26.764 i 


= 0.001 


26.7662 h 


: 0.0004 


26.8 


±0.4 


26.7641 i 


: 0.0001 


B0809+74 


static 


5.735 ± 


0.006 


5.750 ± 


0.007 


10.2 


±0.2 


5.752 ± 


0.001 




model 


5.735 + 


0.005 


5.71 ± 


0.01 


5.5 J 


:0.5° 


5.733 ± 


0.001 


B1133+16 


static 


4.8456 ± 


0.0003 


4.8412 ± 


0.0005 


5.04 


±0.1 


4.8459 ± 


0.0001 




model 


4.8450 ± 


0.0003 


4.8446 ± 


0.0005 


5.15 


±0.1 


4.8451 ± 


0.0001 


B1919+21 


static 


12.4370 i 


: 0.0001 


12.447 4 


: 0.001 


12.8 


±0.1 


12.4373 4 


: 0.0001 




model 


12.4371 d 


: 0.0003 


12.443 4 


: 0.002 


12.5 


±0.1 


12.4370 4 


: 0.0001 



Notes. (fl) DM value computed with step in residuals removed. 



ratio of the observation, but typically we used 4 subbands per 
telescope, corresponding to ~ 12 MHz in the LOFAR data and 
~ 100 MHz in the Lovell data. Each subband was then collapsed 
in frequency and folded in time to create a pulse profile, which, 
unlike the previous analysis, we fitted with Gaussian^ 

The width, height and positions of the Gaussians were free 
to vary, and the best fit values for each of them were recorded 
and used to create a continuous model of how the pulse profile 
evolves as a function of frequency. The fiducial point was chosen 
as either the centre of the main pulse or the midpoint of two 
brightest c omponents d epending on the morphology of the pulse 
profile (see lCralll970h . 

This process is shown for PSR B1919+21 in Figure [3] We 
used power laws to model the frequency-evolution of each of the 
fitted parameters, as they are the simplest functions which fit the 
data well. However, as can be seen in Figure [3]the fits were not 
perfect. The separation of the peaks follows a power law well (al- 
though the exponent is positive, contrary to what is expected by 
radius-to-frequency mapping, see Section 17.2.4) . but the widths 
of the components are not well described by power laws. This is 
due to the overlap between components. In the overlapping re- 
gion it is unclear how much power belongs to each component, 
and the solution which is achieved from the fit is not unique. The 
ratio of the peaks also could not be fitted with a single power law, 
which was the case for all of the pulsars in our sample. Instead 
we used a power law for each frequency band. This is because 
both of the components have different spectral profiles (which 
may contain one or more breaks) and as we don't have abso- 
lute flux values, it is unclear whether one component is getting 
brighter or the other is getting fainter. Scintillation across the dif- 
ferent frequency bands could also contribute slightly to the dis- 
continuities in the component amplitudes, although the effects of 
scintillation will be small in the average pulse profiles because 
of the long integration times used in these observations. 

Our model of the pulsar was then used to create a template 
for each frequency channel, which was cross-correlated with the 
data to get TOAs using pat in the same way as the single tem- 
plates were. 

The right hand side of Figure [T]shows the residuals obtained 
from using the templates derived from this frequency-dependent 
model. Using these model-based templates for timing improves 
the residuals and reduces the systematic errors which were seen 
in the initial residuals. With the frequency-dependent templates 
the residuals show far less systematic trends and agree to within 



6 As noted earlier, for these pulsars Gaussian functions are almost 
identical to von Mises functions, so the different functions used for fit- 
ting the static and model-based templates has no effect on the timing 
residuals. 



the error bars on a single DM in most of our observations (see 
Figure |4] and Table 0. PSR Bl 1334-16 seems to show some 
structure in its residuals. Its DM at ~ 1400 MHz is more than 
3 sigma from its DM at low frequencies and this cannot be re- 
moved by using a frequency-dependent template. This is possi- 
bly due to the finite sampling of the data and template which 
means that at the highest frequencies the pulse profile of PSR 
Bl 1334-16 (whose components are the narrowest of all the pul- 
sars in our sample) is barely resolved. Subsequent observations 
taken at Jodrell Bank, with higher time resolution, show no such 
structure. None of the other pulsars in our sample were affected 
by this problem. 

We will return to the issue of the component evolution in 
pulsars B03294-54, B08094-74, Bl 1334-16 and B1919+21 and 
its impact on pulsar timing (Section|7J, but first we will consider 
the influence of the ISM. 



4. Simulations of Non-dispersive Delays 

Before we started looking for extra structure, the data had al- 
ready been fitted for DM and a phase offset between LOFAR 
and the Lovell telescope. The size of the error bars from the 
cross-correlation are also very different in each of our observ- 
ing bands. For these reasons, it was necessary to perform sim- 
ulations to determine how sensitive the data truly are to these 
frequency-dependent effects. 

To simulate the effects of the ISM, we added a delay to 
the folded profile in each subband of our data according to a 
v~ 4 scaling law (most second-order ISM delays scale roughly 
as v~ 4 , see Section 01. This data was then cross-correlated with 
a frequency-dependent template, fitted for dispersion measure 
and a jump between the LOFAR data and the L-band data, and 
displayed in the same way as the real data to produce residuals 
which had an artificial delay added to them. 

We used the least squares method to fit the residuals with a 
y~ 4 power law. The chi-squared value of the fit and the number 
of degrees of freedom were computed, and used to compare the 
likelihood that the structure in the residuals was caused by a v~ 4 
delay with the likelihood that it was caused by chancel We re- 
duced the magnitude of the delay which was introduced to the 
data until the residuals had an equal probability of being caused 
by chance as being caused by an ISM-like power law. This delay 
was set as the upper limit on the magnitude of this effect in our 
timing residuals. This process is demonstrated in the left panel 
of Figure|5] From these simulations, we determined that our sen- 
sitivity to steep, negative, frequency-dependent power laws is 

7 This was done numerically, using the chi-squared calculator, see 
http://www.fourmilab.ch/rpkp/experiments/analysis/chiCalc.html 
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Table 4. Upper limits derived from simulations. 



PSR 


v 4 delay (ms) 


v u b delay (ms) 




at 48 MHz 


at 180 MHz 


B0329+54 


1.95 


0.65 


B0809+74 


3.84 


1.28 


B1133+16 


1.05 


0.35 


B1919+21 


0.84 


0.28 



dominated by the error bars associated with the LBA observa- 
tions. 

In the same way, we also simulated th e effects of an 
aberration/retardation-like delay. ICordesI (119781) showed that the 
emission height, r(v) oc y~ 2jr , where x is a constant which de- 
pends on the radius-to-frequency mapping of the particular pul- 
sar. The author then used data to show that typically x falls in the 
range 0.21 -0.55. 

The delay in the pulse arrival time is proportional to the 
height at which it is emitted (see Equation IT9b. From this, we 
can deduce that Af should scale with v ' 4-1 . For our simulations, 
we elected to use the average value of x found by Cordes, so 
chose a power law which scales with v° 6 . 

Again, we compared the likelihood of a fitted function and 
no function, and reduced the magnitude of the effect until it 
was undetectable in our timing residuals (see the right panel of 
Figure [5}. Remarkably, the fitted jump has very little effect on 
our ability to detect aberration and retardation effects, and be- 
cause the error bars in our high frequency data are much smaller 
than those at low frequencies we are, in fact, much more sensi- 
tive to high frequency delays. If the power law is steeper than 
y o.6 our sens itivity to these effects increases. Table |4] shows the 
determined upper limits on high and low frequency delays in our 
data. 



5. ISM Effects 

5.1. Impact on Pulsar Timing 

There are a number of second-order effects caused by the ISM 
which could induce additional delays in our data. These effects 
scale strongly with wavelength and are at their strongest at low 
frequencies. As shown by our simulations, the lack of any re- 
maining structure in the residuals indicative of unmodelled ef- 
fects shows there is no significant indication of super-dispersion, 
refraction, anomalous dispersion or frequency-dependent DM 
variations. We can, however, still use these observations to place 
important limits on the magnitude of some of these effects, at 
least for these four sources. 

Taking the maximum unmodelled (i.e. non-dispersive) ISM 
delay in the LBAs and extrapolating it to higher frequencies 
gives an indication of the impact of the ISM o n pulsar timing 
projects like Pulsar Ti ming Arrays (PTAs, see [Roman il l 19891 
Foster & Backer! IT990I) . Figure [6] shows the magnitude of ISM 
effects as a function of frequency for two pulsars with the largest 
(B0809+74) and smallest (B1919+21) deviations from white 
noise in our sample. The largest possible delay in our data (see 
Table |4|, scaled with the gentlest gradient (i.e. the largest pos- 
sible deviation) still corresponds to only ~ 50 ns at 1400 MHz. 
Although this figure will change significantly for every line-of- 
sight, depending on the parameters of the ISM (and potentially 
when the observation was taken) this upper limit is only roughly 
half o f the ~ 100 ns pre cision required for the first generation of 
PTAs dJenet et al.ll2005l) . 



ICordes & Shannon! (1201 Oi and references therein) describe a 
number of frequency-dependent delays, which are caused by the 
distribution and number of free electrons along the line-of-sight. 
Combining these relations with our timing residuals allows us to 
constrain some properties of the ISM (see Table [5}- 

5.2. Frequency-dependent DM Variations 

Frequency-dependent DM variations are caused by scattering 
by the ISM. The diameter of the scattering disk (the region in 
the sky where scattered radiation is received from) increases as 
~ v -2 ' 2 . The DM, which is an average over the scattering disk 
weighted by the flux received from a given direction, therefore 
also varies as a function of frequency (Figure|7). 

If we assume that all of the scattering material is concen- 
trated in a thin screen ha lf-way along the line-of- sight, the vari- 
ation in DM is given by (Cordes & Shannon 201(J: 

6DM = 2.82 x 10 4 D' 5/ V I1/6 SM pc cm" 3 , (6) 

where D' is the distance between the pulsar and the scatter- 
ing material (in kpc), v is the observing frequency in MHz and 
SM is the Scattering Measure along the line-of-sight in units of 
kpc irT 20/3 . This corresponds to an RMS variation in time delay 
of: 

AfDM.v = 3.79 x 10 4 D' 5/ V 23/6 SM s. (7) 

This allows us to constrain the SM along the line-of-sight (see 
Column 4 of Tabled. The SM must be fairly large for this effect 
to be detectable and all of the pulsars which we observed are at 
low DM (< 30 pc cirT 3 ) and show very little scattering (with the 
exception of PSR B0329+54, no scattering is visible in any of 
the pulse profiles). Our upper limits are about 2- 3 orders of mag- 
nitude above those from the NE2001 model of ICordes & Laziol 
(120021) . 

To detect this effect in PSR B0329+54 (assuming the SM 
from NE2001) would require a timing residual RMS of 10 yus, so 




Frequency (MHz) 

Fig. 6. These curves show the upper limits on the size of de- 
lays from second-order ISM effects extrapolated up to higher 
frequencies. The dark grey area is for PSR B0809+74, which 
had the largest RMS residuals in our timing fits; the light grey 
area is for PSR B1919+21, which had the smallest RMS resid- 
uals. They are scaled with v~ 3 and v~ 4 which are the lower and 
upper bounds for scaling of the second-order ISM effects. 
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SCATTERING SCREEN 

Fig. 7. The size of the scattering cone is larger at low frequencies 
(dark line) than at high frequencies (dashed line). The DM is a 
flux-weighted average over the area of the scattering disk so it 
varies with frequency as the scattering cone changes size. 

it is not surprising that this effect was not detectable in our data. 
However, it is possible that this effect could be detected in the 
future with LOFAR using pulsars with a higher SM (for example 
PSR J2044+4614 or PSR B2036+53, both of which are closer to 
the Galactic plane). The LOFAR data presented here is far from 
full sensitivity, and, when the data are fully calibrated, and all 24 
of the stations in the LOFAR core can be combined coherently, 
the telescope's sensitivity will be improved by at least an order 
of magnitude. 

5.3. Refraction 

Refraction can also introduce a delay into pulsar timing data at 
low frequencies. If there is a gradient in the interstellar electron 
density perpendicular to the line-of-sight, rays of light are re- 
fracted, and bent through so me angle 6 r . For a thin screen this is 
given bv lClegg et all (119981) : 



A 2 , 



2jt dx 



N e (x), 



(8) 



where r e is the classical electron radius and N e is the column 
density of electrons along the line-of-sight. At lower frequen- 
cies, light is bent through a larger angle. The larger the angle, 
the further the light must travel and the longer it takes to arrive 
at t he observer (see F igure [8). 

ICorde s & Shannon (2010) showed the delay between a 
straight path and a refracted path is given by: 



1 i 

Af re f = —D eB 6 2 . 

2c 

where D e ff is a characteristic distance, given by: 



Deff = 



D-D s 
DJD 



(9) 



(10) 



D is the distance between the observer and the pulsar and D s 
is the distance between the observer and the scattering screen. 
Assuming the line-of-sight is dominated by one screen (which is 
a re asonable assumptio n given the evidence in IStineb ring 2006 
and|Brisken et al. 2010, and that our sample is quite nearby), the 
thin screen model is valid. By substituting 9 r from Equation [8] 
into Equation [9] and rearranging for v we can derive an expres- 
sion for the time delay between refracted and unrefracted rays: 



Af ref = 3.47 x 10 1 



jDeS ( d 

v 4 \dx 



DM 



Fig. 8. A gradient in the electron density perpendicular to the 
line-of-sight causes rays to be bent. The size of the angle which 
the rays are bent through depends on the frequency of the radia- 
tion. Low frequencies (dark line) are bent through a larger angle 
and so must travel along a longer path to reach the observer, de- 
laying them with respect to higher frequencies (dashed line). 



where ^DM is the average gradient in the DM perpendicular 
to the line-of-sight in pc crrT 3 AU -1 , v is in MHz and D er r is in 
kpc . This number can be constrained from the lack of structure 
in our timing residuals and can be used as an indicator of how 
much the DM of a pulsar is likely to change. 

It follows from the derivation above, that the -r-DM value ob- 
tained holds over distances dx of the size of the scattering cone. 
At these frequencies the scattering cone for a typical (nearby) 
pulsar is on the order of a few AU (for B0329+54 the scatter- 
ing cone is 10 AU at half the distance to the pulsar), so the most 
natural units to use for this quantity are pc cirT 3 AU~' . The unit 
pc cirT 3 ALT 1 also corresponds to how much the DM of a pul- 
sar at a distance D travelling 250 km s _1 will change in approxi- 
mately one week. It should be noted however, that the number is 
an average over the entire scattering disk, and so is not sensitive 
to small scale anisotropies. Upper limits are shown for j^DM for 
the pulsars in our sample in Column 5 of Table [5] 

The main source of error in this number arises from D e ff as 
the distance to the scattering screen is unknown. The function 
is very sensitive to nearby scattering screens (D s < 0.25D) but 
is not very sensitive to distant scattering screens {D s > 0.75D). 
It is, however, approximately constant for 0.25 < D s < 0.75, 
varying by only a factor of ~ 3 in this range. In our calculations, 
we assume that the scattering screen is roughly half way along 
the line-of-sight and D er i « D. 

5.4. Anomalous Dispersion 

By modifying Equation Q] to include gyro-motion in a mag- 
netic field and electron collisions, and including the quartic term 
in the Taylor expansion of l/v g , we can write a more general 
version of the plasma dispersion law, which can be expressed 
as the sum of three functi ons (Equations Qj] [15] and Q~6] see 
iPhillips & Wolszczarifl992l) : 



Af DM = At y , + Ar v , + At v 



The j\ term: 



(ID Aty 



DM 



2.41 x 10-V 



(l+#herm)s 



(12) 



(13) 
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is the normal dispersive delay term multiplied by an additional 
factor which depends on the ratio of the thermal velocity of the 
electrons to the speed of light, /Jtherm : 



Aherm — \ T • 



kT 



(14) 



Where k is the Boltzmann constant, T is the temperature of the 
plasma and m is the mass of an electron. This means that if the 
temperature is non-zero it adds a small contribution to the disper- 
sive delay. In practice this is a very small number and is absorbed 
into the v~ 2 law making the DM slightly higher than it would be 
normally, but not altering the power law of the delay. Even an 
unrealistic change in temperature as large as 100 K only modi- 
fies the apparent DM by 0.01%, and would be indistinguishable 
from the normal dispersive delay. 

The y2 term is because of the plasma being weakly magne- 
tized. For a circularly polarised wave: 



Aty 



0.0286RM 

73 



(15) 



where RM = JT n e B\\Al is the Rotation Measure in rad m 2 . This 
only affects circularly polarised sources and is generally negligi- 
ble. The RM along a given line-of-sight is closely related to the 
SM and the DM, so when a pulsar's RM is large enough to make 
this effect detectable, scatter broadening makes the pulsar unde- 
tectable. An RM of 280 would be needed to produce a detectable 
delay of 1 ms at 20 MHz. The scattering time for a pulsar with 
this RM would be ~ 4 s, which, in most cases, would render 
the pulsar undetectable (assuming RM ~ 5 X DM, the average 
from the ATNF pulsar catalogue Manchester et al. 2005, and the 
empirical scattering law from lBhat et al.ll2004 ). 

The 73 term depends on the Emission Measure, EM = 

rD -) 

J rigdl along the line-of-sight: 



Aty 



EM 
4-Ov 4 



(16) 



where EM is given in pc cm 6 . For a uniform distribution of 
plasma between the pulsar and observer the contribution from 
the v~ 4 term is small even at low frequencies and is proba- 
bly not detectable. However if the electrons are arranged in 
a clumpier distribut ion along the line-of-sight (as discussed in 
iKuz'min et ai] [2008) this term becomes larger. 

Whilst both the j3 t i ierm and 72 terms are too small to detect, 
the fact that no delay is observed in our data can be used to con- 
strain the EM along the line-of-sight through the 73 term. The 
upper limits are given in Column 6 of Table [5] Although they are 
currently around five orders of magnitude hig her than the values 
predicted by NE2001 dCordes & L azio 2002), using this method 
on observations at lower frequencies and with higher sensitiv- 
ity, could provide a new way to measure the EM along a given 
line-of-sight. 



6. Magnetospheric Effects 

In addition to the delays caused by the ISM, there are potential 
sources of delay from within the pulsar magnetosphere itself. 
We can use the absence of any additional delays in our timing 
fits to constrain the composition of the magnetosphere and the 
emission height above the magnetic poles. 



6.1. Aberration and Retardation 

According to radius -to-frequency m appin 

(iRuderman & Sutherland! Il975t iManchester & Tavlorl Il977l 
Cordesl ll978l) . high frequency emission is thought to originate 
close to the neutron star whilst low frequency emission comes 
from higher in the pulsar magnetosphere. The range of emission 
heights at different frequencies can be obtained, independently 
of the radius-to-frequency mapping mo del, through t he effects 
of differential aberration and retardation (1Cordeslll978b . 

Retardation is the time delay caused by the difference in 
path length from the different emission sites. For emission which 
originates at an altitude r m ; n < r max , the time delay between the 
two paths is given by: 



At n 



(17) 



Aberration is caused by the co-rotation of the magnetosphere, 
which bends the radiation beam towards the direction of rotation 
and therefore causes pulses to arrive earlier than they would if 
they were to travel along a straight path. Aberration increases at 
larger radii, so the result is to delay pulses from low altitudes by: 



At 



ab 



sin a At,-, 



(18) 



where a is the angle between the pulsar's magnetic and rota- 
tional axes. Combining the two effects gives the total time delay 
for a given radius: 



At 



A/R 



(1 + sin a)- 



(19) 



In our data, assuming the fiducial points in our frequency- 
dependent models are aligned correctly, there is no remaining 
structure in the residuals that has not been successfully mod- 
elled by a quadratic frequency dependence. This can be used to 
constrain AR = r max - r m ; n , the distance between the heights at 
which different frequencies are emitted (see Table |6). 

For pulsars B0329+54 and Bl 133 +16, which ex hibit typi- 
cal 'conal' behaviour (see for example Rankin 1983a), it is also 
possible to use radius-to-frequency mapping to determine upper 
limits on the absolute height of emission from the pulsar sur- 
face. At low frequencies radius-to-frequency mapping suggests 
that pulse profiles are broadened as the star's dipolar magnetic 
field lines move further apart high in the magnetosphere. For a 
neutron star with a dipolar magnetic field the ratio of the widths 
of the profiles (8 V 2 > 8 V \) at two frequencies ( vi > yg) can be 
related to the emission heights by the equation dCordeslll978l) : 



AR 

'"min 



— I -1. 



(20) 



This equation can be used in conjunction with the values of AR 
derived earlier to determine upper limits on r max , the absolute 
height of the 40 MHz emission. This analysis was not performed 
for pulsars B0809+74 and B1919+21, as our observations of 
their pulse profile evolution do not agree with the standard pic- 
ture of radius-to-frequency mapping (see Section |7). 

Our limits agree well with previous papers, such as those by 
ICordesl(fT978h : lMatese & Whitmirel(ll980l) : lKaruppusamv et al.1 
(1201 ll) . who also failed to detect the effects of aberration and re- 
tardation in low frequency pulsar timi ng data. It is a l so int erest- 
ing to compare our findings to those of Kra mer et al.1 (1 1997 b who 
performed a similar analysis at high frequencies (1.4-32 GHz) 
and found that r max < 310 km for PSR Bl 133+16 and r max < 
320 km for PSR B0329+54. 
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Table 5. Derived values of the dispersion measure (DM), and upper limits on the scattering measure (SM), perpendicular gradient 
in dispersion measure (^DM) and emission measure (EM) derived from the analysis described in the text. The values of DM are 
shown here only for comparative purposes, and correspond to the average DM derived earlier from the frequency-dependent model. 
For more information, see Table [3] 





Delay 


DM 


SM 


£dm 


EM 




(ms) 


(pc cnT 3 ) 


(kpc irr 20/3 ) 


(pc cnr 3 AIT 1 ) 


(pc cm" 6 ) 


B0329+54 


1.95 


26.764 


< 0.25 


< 5.3 x 10" 


< 42000 


B0809+74 


3.84 


5.733 


< 1.02 


< 1.2 X 10~ 5 


< 82000 


B1133+16 


1.05 


4.845 


< 0.32 


< 6.7 x 10" 6 


< 22000 


B1919+21 


0.84 


12.437 


< 0.155 


< 4.4 x 10" 6 


< 18000 



Table 6. Constraints on emission heights for frequencies be- 
tween ~ 40 and ~ 180 MHz from aberration/retardation argu - 
ments. The a values are taken from lLvne & Manchester! dl 988). 
r max was not calculated for pulsars B0809+74 and B1919+21 
as the frequency evolution of their pulse profiles is not consis- 
tent with a dipolar magnetic field and simple radius-to-frequency 
mapping. 





A?A/R 


a" 


AR 


'"max 




(ms) 


o 


(km) 


(km) 


B0329+54 


< 0.65 


30.8 


< 128 


< 183" 


B0809+74 


< 1.28 


0.0 r 


< 384 




B1133+16 


< 0.35 


51.3 


< 59 


< 110 


B1919+21 


< 0.28 


45.4 


< 49 





Notes. (n) Errors were not given for a in lLvne ^ Manchester ( 198j|), so 
errors on the emission heights cannot be derived. See text for a discus- 
sion. (A) Emission height of inner cone. (c) 0.0 was used to provide an 
upper limit because alpha is unavailable for B0809+74. 



Our limits of r max < 110 km for PSR Bl 133+16, 
significantly improve upon the previous low freque n cy (< 
200 MHz) limits for P S R Bl 133+16 set by ICordesI (119781) 
and iMatese & Whitmirel d 19801) . who found r max < 630 km; 
and lKaruppusamv et al.l (1201 II) who improved upon this, finding 
''max < 560 km. This is predominantly because we have a large 
fractional bandwidth and high sensitivity at low frequencies. Our 
limits are also improved by the frequency-dependent models that 
were used, which allowed us to test how well our fiducial points 
fit the data, reducing the uncertainty in AfA/R significantly. 

The uncertainties in our measurements are dominated by 
the uncertainties associa ted with a, which unfortu nately, are not 
well constrained (see lEverett & Weisbergl [200 1[ for a discus- 
sion). No uncertainties on a were given in iLvne & Manchester! 
( 1988), and because of this, the uncertainties on AR and r max are 
impossible to determine definitively. However, as a only appears 
in Equation[19]through a factor of (1 + sin a) our measurements 
should be within a factor of two of the true value. 

The implication of these limits is that pulsar emission from 
all radio frequencies is produced inside a very small region of 
the magnetospher^- All of the radio emission from Bl 133 + 16 
comes from within 1 1 stellar r adii (using the can onical neutron 
star radius of 10 km, as used in lKramer et al.ll997l) . a fact which 
could have implications for future models of the pulsar magne- 
to sphere. 



8 Our calculation of absolute emission height (r max ) assumes that the 
emission comes from dipolar magnetic field lines emanating from the 
polar cap, although the calculations of the range of heights, ( AR) is valid 
for any geometry. 



6.2. Magnetospheric Propagation Effects 

The pulsar magnetosphere is a complex system with strong mag- 
netic fields and high concentrations of relativistic charged parti- 
cles. As high frequency emission is thought to originate closer 
to the neutron star surface, it has more of the magnetosphere to 
travel through. This means that, under the assumption of radius- 
to-frequency mapping, we might expect to measure a slightly 
higher value for DM at high frequencies than at low frequencies, 
changing the way dispersion delay scales with frequency. 

We see no evidence for any deviation from a v~ 2 power law 
in our data, suggesting that dispersion from within the magneto- 
sphere is either not present, too small to detect, or indistinguish- 
able from the cold plasma dispersion law (at least in the fiducial 
components). This suggests that either the column density of the 
plasma in the magnetosphere is too small to cause refraction, or 
emission can somehow escape the magnetosphere without be- 
ing refracted. This could be because the emission of the fiducial 
component propaga te via the extraordinary mode, which does 
not suffer refraction (Bar nard & Aronsll986b . In the pulsar mag- 
netosphere, the electrons are very tightly bound by the magnetic 
field lines and cannot move transverse to their direction. This 
means that photons cannot excite them, effectively making them 
invisible and setting their refractive index to 1 (see for example 
Michellll991h . It should be noted that other modes of propaga- 
tion, which can be refracted (ordinary modes) are also possible 
in the magnetosphere, although from this evidence, they do not 
seem to be present in the fiducial components. 

6.3. Super-dispersion 

Super-dispersion was proposed by Shitov & Malofeev (19851) af- 
ter they observed that the DMs of low frequency pulsar observa- 
tions seemed to be systematically higher than the DMs of the 
same pulsars at high frequencies. 

T hey explained the delay by magnetic sweepback dShitovl 
[1981 . in which a pulsar's magnetic field lines get bent back- 
wards in the opposite direction to the spin of the neutron star. 
This causes emission from higher up in the magnetosphere (at 
low frequency) to reach us slightly later than the correspond- 
ing emission from closer to the neutron star surface (high fre- 
quency). This m odel was supported by further evidence from 
Kuz'min (1986) who observed super-dispersion in eight more 
pulsars and observed a correlation betw een 1 IP and th e observed 
super-dispersive delay and also Shit ov et al.l ([1988) who ob- 
served the effect in pulsars B0834+06, Bl 133+16, B1508+55, 
B0832+26 and B1642-03, and also noted that the delay ap- 
peared to be time v ariable in Bl 133 + 16 and B0809+74. 

In a later paper (iKuz'min et alJ2 008) this was cast into doubt 
as the super-dispersion in Crab giant pulses corresponded to 
more than 1 period, a delay which cannot be explained by the 
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twisting of magnetic field lines in the pulsar magnetosphere. 
Super-disp ersion was also no t seen by iPhillips & Wolszczanl 
( 1992{) and lAhuia et alj d2005l) . who did, however see an excess 
delay at high frequencies. 

We see no evidence for extra delays at low frequencies and 
can place a limit on the super-dispersive delay of < 1 ms at 
40 MHz for the pulsars in our sample. We speculate that the 
delay which was seen in PSR B0809+74 was actually due to 
either pulse profile evolution or an incorrect fiducial point (see 
Section|7]). 

7. Pulse Profile Evolution 

7.1. Impact on Pulsar Timing 

We have shown that it is possible to find an analytical fiducial 
point for each of the pulsars in our sample, which is valid be- 
tween 48 MHz and 1780 MHz, and which satisfies the v~ 2 fre- 
quency dependence expected from the dispersive delay. Building 
a frequency-dependent template around this fiducial point to in- 
clude the effects of pulse profile evolution across the band sig- 
nificantly reduces the systematic errors caused by pulse profile 
evolution, and improves the precision of timing observations. 

Pulse profile evolution can cause systematic errors in pulse 
arrival times which can be on the order of milliseconds for nor- 
mal (slow) pulsars. The size of this effect is largest when th e 
profile evolution is asymmetric, as noted in [Ahuia etahl (120071) . 
but it still plays a role in relatively symmetric profiles like that 
of PSR B0329+54. In reality there are very few pulsars that have 
truly symmetric pulse profile evolution, so in order to obtain ever 
higher timing precision (on the order of microseconds for normal 
pulsars or tens of nanoseconds for millisecond pulsars) it is cru- 
cial to account for the evolution of pulse shape across the band. 

Although the frequency evolution across the relatively small 
bandwidths used up until recently will limit the effect in an indi- 
vidual band, it will still manifest when one tries to combine data 
from more than two frequencies as it cannot be absorbed into a 
fit for a dispersion delay. The problem becomes more acute as 
we search for greater sensitivity by using wider and wider in- 
stantaneous bandwidths. Here the determination of the time of 
arrival itself is affected by the evolution of the profile and the 
dispersion delay. The method presented here provides a way in 
which one can use very wide band data to model the profile suf- 
ficiently to build a template which incorporates all these effects, 
although it remains to find the optimal way to extract a time of 
arrival from these data. 

7.2. Models 

To address the profile evolution induced errors in the residuals it 
was necessary to make frequency-dependent models for each of 
the pulsars in our sample. We were able to model the profile evo- 
lution of the four pulsars over seven octaves of frequency using 
analytic models of Gaussian fits to the data. Where necessary (in 
PSR B0329+54) the models were also made to include the ef- 
fects of interstellar scattering. This was modelled by convolving 
an exponential tail (whose length was fit to the particular fre- 
quency band) with the Gaussian components. Although the mod- 
els used are simple, they describe the shape of the pulse profile 
very well (as shown below), and are very effective in reducing 
the systematic errors seen in our timing residuals as a function 
of frequency. Our timing residuals were used to test the validity 
of each of our models by determining how well the frequency- 
dependent templates remove the different systematic errors at- 



tributed to profile evolution in each observation. The power-law- 
dependencies of all of the fitted parameters in the models are 
given in Table Q 

In all cases, the evolution of the relative amplitudes of the 
peaks was very difficult to fit. In general the peak heights within 
an observing band could be approximated well by a power law, 
but the parameters of each power law varied significantly from 
one observing band to the next, leading to discontinuities in our 
model. This could be because the pulse components do not have 
the same spectral index across the entire frequency range (i.e. 
they have one or more spectral breaks), which could cause the 
observed discontinuities in the 'ratio of peaks' parameter. The 
only means to test this hypothesis requires accurate flux mea- 
surements for the various components and since absolute flux 
calibration is not yet available to us, the present data set cannot 
be used to provide a conclusive answer to this particular ques- 
tion. 

7.2.1. PSR B0329+54 

iGangadhara & Gupta! d2001l) (hereafter GG) model PSR 
B0329+54 as a central core component surrounded by four 
nested conal rings. Our model agrees well with this; we use a 
central core and two cones (five components) because rings 2 
and 4 of the GG model are too weak for us to detect with any of 
our instruments. Figure |9]shows our model. 

The core component (component 3), is the fiducial point for 
the observation and so we have set it to be constant in height and 
position, allowing the other components to vary. The outer cone 
(components 1 and 5) fades at low frequencies, and is too weak 
to see in the LOFAR observations. Interestingly, the two sides 
of the cone fade at different rates, component 5 has a steep spec- 
trum, and is not detected in the HBA or LBA observations (it was 
detected by GG at 606 MHz and 325 MHz), whilst component 1 
fades more gradually, and does not disappear until the LBA ob- 
servations. The inner cone (components 2 and 4) is brighter and 
is visible all the way down to the LBAs (although the scattering 
tail makes it more difficult to model the components here). Again 
there is evidence for the two sides of the cone showing different 
spectral indices. Component 2 fades more than component 4 in 
the high frequency observations, but component 4 fades more in 
the LBA band. 

The widths of all of the components seem to remain re- 
markably constant over the entire frequency range of our ob- 
servations. In fact, PSR B0329+54 can be modelled with con- 
stant component widths from frequencies between 40 MHz and 
8 GHz. It is difficult to say conclusively whether this model re- 
flects a feature which is intrinsic to the pulsar, as the components 
all overlap, making them difficult to model. However a model 
with fixed component widths is simpler and this makes it much 
easier to track how the pulse profile evolves at different frequen- 
cies, how the components move and how their brightnesses vary 
in relation to each other. It is also worth not ing that similar be - 
haviour has been found in the Vela pulsar bv lKeith et alJ (1201 ll) . 

Both cones are very asymmetrical in terms of the relative 
brightness of their two peaks and their relative positions in re- 
lation to the central component (see Figure ITOt. Compared to a 
model where the cones are symmetric around a central compo- 
nent, the outer cone is skewed by approximately 5 degrees to- 
wards earlier pulse longitudes. The leading component moves 
away from the central component with decreasing frequency, 
whilst the trailing component seems to move slightly closer. The 
inner cone is also skewed by approximately 5 degrees, but in 
the opposite direction. The components both move out from the 
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Fig. 9. The model of PSR B0329+54 used for the dynamic tem- 
plate. The pulsar is modelled using five Gaussian components 
(plotted with grey lines). At low frequencies the model is con- 
volved with an exponential function to account for the effects of 
scatter broadening. The final model (including the scatter broad- 
ening) is plotted in black. 




Frequency (MHz) 

Fig. 10. A fit to the relative positions of the components of PSR 
B0329+54, which show a lot of asymmetry. The outer cone (bold 
line) is skewed towards earlier pulse longitudes and both of its 
components fade at very different rates. The inner cone (dashed 
line) is skewed in the opposite direction and again shows differ- 
ent spectral indices for each side of the cone. Models are less 
reliable in the LBAs (< 100 MHz) where the pulse profile is 
affected by scattering. 



main pulse, but component 4 moves out quite quickly, and re- 
places component 5 as the 'postcursor' in the HBAs, whilst com- 
ponent 2 remains in roughly the same place until the top of the 
LBAs, when it begins to move out. 

Our model agrees well with the GG model of a central 
core component surrounded by cones. The cones, however, 
show asymmetrical behaviour which is not expected for a dipo- 
lar magnetic field, the radio emission should come from con- 
centric circles centre d on the magnetic pole (see for example 
lOster & SiebeilT976l> . 

7.2.2. PSR B0809+74 

PSR B0809+74 is a rather controversial pulsar, which has 
shown evidence of an 'absorption feature' (Battel 1981; Rankin 
1983b) and was also one o f the candidates for 'super-dispersion' 
(IShitov & Malofeevlll985L see S ectionlQ). 

Absorption was proposed by B attel et aT1(ll98ll) who noticed 
that the DM found in 102 MHz observations was significantly 
different from the value found from observations at 1720 MHz. 
The reason for this was that their fiducial point was on the lead- 
ing edge of the low frequency pulse profile, and the trailing edge 
of the high frequency pulse profile. It appeared that the low 
frequency pulse profile was 'missing' radiation from the lead- 
ing edge, which they suggested, was removed by cyclotron ab- 
sorption. Further evi dence for thi s model was also provided in 
a subsequent paper (Battel 1981), which found that th e profile 
of B080 9+74 gets significantly narrower below 1 GHz. Rankin 
( 1983b) found similar absorption features in at least eight other 
pulsars. 

PSR B0809+74 has two overlapping components, which are 
normally thought to be conal. In accordance with this thinking, 
we fitted the data from the simultaneous observations (marked 
with arrows in Figure [PD with two Gaussian components and set 
our fiducial point as the midpoint of the profile. This model pro- 
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Fig. 11. The model used to produce the dynamic template of PSR 
B0809+74. The model consists of two Gaussian components. 
The peak of the narrower component is the fiducial point of the 
observation and the broad component drifts through the pulse 
profile. The two components and the final model are plotted in 
black, and compared to data, which is plotted in grey. The si- 
multaneous observations (used to create the model) are indicated 
by aiTows. Pulse profiles at 410 MHz, 606 MHz and 925 MHz 
are from the EPN database and the low frequency (10-60 MHz) 
pulse profiles are from a recent observation with the LOFAR su- 
perterp. 
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Fig. 12. The position of component 1 relative to component 
2. The data follow a single smooth power law (which scales 
as v -°- 43±0 06 ) suggesting that the component drifts through the 
pulse profile. 

duced large systematic errors in the TOAs at different frequen- 
cies. On a closer examination of the profiles, the reason for the 
timing errors became apparent, the separation of the components 
cannot be modelled as a simple power law. In the observations 
at 1400 MHz the components get closer together as frequency 
decreases, whereas in the low frequency data they move further 
apart. 

In a second model, we used three components, one central 
component, a precursor and a postcursor. The narrower, central 
component is taken as the fiducial point of the profile. At high 
frequencies the precursor (component 1 in Figure [TT] moves to- 
wards the central component. Somewhere in the frequency range 
200-1000 MHz (which was not present in the simultaneous ob- 
servations) this component fades. Then, in the low frequency 
data, the third component appears and begins to move away from 
the central component, towards later pulse phase. 

The precursor and the postcursor in this model have the same 
width and their positions can both be modelled by a single power 
law (see FigurefT2l. This suggests that the two components may 
instead be a single component, which drifts through the pulse 
profile. 

In our final model (see Figure [TTJ, we used two Gaussian 
components, a narrow component (component 2), which is the 
fiducial point of the pulse profile, and a broader component 
(component 1), which starts as a precursor in the high frequency 
observations and drifts through the pulse profile, arriving at a 
later phase (as a postcursor) at low frequencies. Using the nar- 
row component as a fiducial point removes the systematic errors 
from TOAs, which provides strong verification of this model. 

Further evidence in favour of this model comes from 
archival pulse profiles from the E uropean Pulsar Network 
(EPN) database dGould & Lvne1[T99 8fl These pulse profiles (at 
410 MHz, 606 MHz and 925 MHz) are also shown in Figure [TT1 
along with an interpolation of our model to these frequencies. 
We have allowed the relative heights of the components to vary, 
but their positions and widths are determined by our model. 
Without prior knowledge about this frequency range, our model 
accurately predicts the shape of the pulse profile. 

9 http://www.mpifr-bonn.mpg.de/old_mpifr/div/pulsar/data/ 
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The fitted solution for the profile at these frequencies (where 
the two components overlap), is not unique. However, it is ob- 
vious from considerations of the timing residuals that the mid- 
point of the two components is not the fiducial point of the pulse 
profile and the fact that such a simple model can explain the 
observed profile evolution of the pulsar over such a broad fre- 
quency range is compelling. 

We also compare our model to some more recent LBA ob- 
servations taken using the LOFAR superterr£3 (LOFAR obser- 
vation ID: L30803). The data quality is significantly improved 
because the superterp has roughly three times the collecting area 
of DE60 1 , and the delays between the dipoles have recently been 
calibrated. The pulse profiles are made from 6 MHz segments of 
bandwidth between 15 MHz and 57 MHz. Our model accurately 
describes the two components down to roughly 40 MHz, where 
the broader component begins to move away from the central 
component more slowly. This could be due to a mode change or 
some more complex pulse profile evolution at low frequencies 
(perhaps betraying one of the magnetospheric effects discussed 
above). 

The polarisation of the EPN profiles also shows some ev- 
idence that one of the components in our model is linearly po- 
larised (see FigureQj). The first component in the data is linearly 
polarised at high frequencies, but as frequency decreases the lin- 
ear polarisation moves towards later pulse longitudes. The early 
LOFAR polarisation data (obs ID: L24117), shown in the bot- 
tom panel of the figure, shows that the polarisation has moved 
towards the latter portion of the pulse profile at 136 MHz, and 
arrives after the main pulse. 

One argument against this model is the step in subpulse 
phase which is observed at 1380 MHz (see for example 
lEdwards & Stappersl2003al) . If the precursor at 1380 MHz is the 
same component that appears on the trailing edge of the profile 
at low frequencies, it is difficult to explain what happens to the 
phase jump, which is not observed at 328 MHz. The mechanism 
by which the two components are seen to move through each 
other, is also still unknown. It could be due to a retardation effect 
which is only present in one component, or refraction from the 
magnetosphere, but is not immediately obvious what could be 
causing the profile to evolve as it appears to. Further investiga- 
tion into the single pulses of PSR B0809+74 at low frequencies 
could help to answer these questions. 

It is interesting to note that the fiducial poin t in ou r model 
matches the fiducial point used by iBartel et ail (Il98ll) . which 
led the authors to speculate that part of the low frequency pulse 
profile was missing. Our model, shows that the radiation is not 
missing, but has been displaced somehow, to later pulse longi- 
tudes. Our model also elegantly explains the narrowing of the 
pulse profi le (also attribute d to absorption) discussed in iBartell 
(Il98lh and lRankinl (Il983bl) . The cumulative pulse starts at high 
frequencies as two fairly distinct components, the components 
get closer together as frequency decreases, reducing the appar- 
ent width of the profile. At around 400 MHz, the profiles are ex- 
actly on top of each other, and the pulse width is at a minimum. 
Below this frequency, the broader pulse continues to move to- 
wards later pulse phase and the profile width begins to increase 
again, reproducing the shape of the absorption. 

Super-dispersion can also be explained by this model, if the 
centre of the two components was used as a fiducial point for 



10 The superterp is a grou p of six core station s, whose signals can be 
combined coherently (see Stap pers et al.ll20lH for more details), cur- 
rently the most sensitive LOFAR observing configuration for pulsars 
and fast transients. 
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Fig. 13. Linear polar isation profiles of PSR B0809+74 between 
925 an d 328 MHz dGould & Lvnel [19981: lEdwards & Stappersl 
I2003bl) . and LOFAR polarisation profile at 136 MHz (black 
lines) plotted along with the Stokes I profiles at each frequency 
(grey lines). The polarised component moves from the position 
of the leading component of the profile towards later pulse lon- 
gitudes, tracing the broad component of the pulse profile. 



the pulse profile, instead of the core component, as in our model 
then the pulse would look like it arrived later than expected at 
low frequencies. 

7.2.3. PSR B1 133+16 

PSR Bl 133+16 is one of the prototypical exa mples of a 'well 
resolved conal double' profile (Rankin Hl983al) . It is one of the 
brightest pulsars in the Northern sky, so it has been widely stud- 
ied and is often used as evid ence in favour o f radius-to-frequency 
mapping (see for example iThorsettl 1 1 99 ll) . The separation be- 
tween its components shows a continuous increase with decreas- 
ing frequency, which is thought to trace the dipolar shape of the 
pulsar's magnetic field. 

This is exactly what we see in our model. It has two strong 
components (components 1 and 3) separated by bridge emission 
(component 2) (see Figure fT4l . The fiducial point is the mid- 
point between components 1 and 3, which has been defined to 
be the peak of the bridge emission in this model. The conal com- 
ponents move further apart at lower frequencies, scaling with a 
power law ~ y~ °- 62 . This is consistent wi th the exponent found in 
IThorsettl (I 1 99 ll) and lXilouris et al.l(ll996l) . who found exponents 
of -0.50 and -0.71 respectively . The exponent is however , sig- 
nificant l y low er than the value in Karuppusamy e t alj d201 ll) and 
ICordesI (1 19781) . who both found a power law ~ v . 
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Fig. 14. The model used to produce the dynamic template of PSR 
Bl 133+16. The source is modelled as three Gaussian compo- 
nents (plotted in grey): two conal components and one which is 
attributed to bridge emission. The final model is also plotted in 
black. 



This is because their power law fits did not include the con- 
stant term, A# m j n , the width of the pulse profi le at the surfac e 
of the neutron star. This term was proposed by iThorsettl (119911) . 
who found that the separation of the components tends towards 
a constant value at high frequencies. If we do not include this 
term in our fits, we find a power law of v~° 2 , whic h is roughly 
consis tent with those of Karuppus amv et aD (1201 ll) and|Cordes 
dl978l). For our analysis, we used the form of the power law from 
Thorsett (1991), which provides a better fit to the data. 

The width of the conal components as a function of fre- 
quency has been a subject of interest in the past. lMitra & RankinI 
(2002) found that the component widths remain constant be- 
tween 40 and 3000 MHz. We also see evidence of this in our 
model above 80 MHz, altho ugh the component wid ths begin 
to broaden below this value. iMitra & RankinI (120021) also no- 
ticed this broadening, and attributed it to dispersive smearing 
across a frequency channel or scattering from the ISM. In our 
observations, the dispersive smearing at 48 MHz (across a sin- 
gle 12 kHz channel) is ~ 1.5°, which is enough to explain the 
observed broadening in the profileEl. However, even disr egard- 
ing this low frequency broadening Mitra & RankinI (120021) found 
that the spacing between the components between 100 MHz and 
10 GHz changes too rapidly to be caused by a dipolar magnetic 
field. 

PSR Bl 133+16 is the pulsar which is most consistent with 
radius-to-frequency mapping of all the pulsars in our sample. 
However, in Section [6] we showed that its emission is confined 
to a very narrow region in the magnetosphere (<59 km) which 
is incompatible with the standard radius-to-frequency model. 
Radius-to-frequency mapping assumes that the emission at a 
given emission height traces the last open field line in the pulsar 
magnetosphere. From ge ometrical arguments (see for example 
ILorimer & Kramerll2005l) it is possible to write the opening an- 
gle of the last open field line of a non-relativistic dipole, p, as a 
function of emission height, r. 



P = 86° , 



(21) 



where R^c is the radius of the light cylinder of the pulsar, p is 
the maximum value that it is possible for components to be sepa- 
rated by as each component originates from one side of the dipo- 
lar field. The maximum increase in component separation pre- 
dicted by this equation for PSRB1133+16 over the height range, 
AR, derived earlier is 2.8°. In the pulse profile of PSR Bl 133+16 
the components are seen to move apart by 6.3±0.5° between 
17 80 MHz and 48 MHz. This, coupled with the further evidence 
bv lMitra & RankinI d2002l) suggests that PSR Bl 133+16 cannot 
be explained through radius-to-frequency mapping in a simple 
dipolar magnetic field. 

7.2.4. PSR B1919+21 

PSR B 19 19+21 w as the first pulsar ever discovered 
dHewish et all 1 19681) and is often referred to as a classic 
radio pulsar. In fact, PSR B1919+21 disagrees with the classic 
picture of a pulsar in almost all aspects of its pulse profile 

ev olution. 

iLvne et al] d!971l) showed that the two components of the 
pulsar get closer together as frequency decreases in the range 
150-3000 MHz, which is the opposite to what one would expect 



11 The half power width of component 1 increases from 1.9° at 
72 MHz to 3.6° at 48 MHz, and component 2 increases from 2.1° at 
72 MHz to 3.4° at 48 MHz. 
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from radius-to-freq uency mapping. This result was confirmed by 
Siebe r et all (1 19751) who also showed th at the pulse profile seems 
to get broader again below 150 MHz. iMitra & Rankinl d2002l) 
found that the width of the profile and separation of the compo- 
nents was approximately constant in their observations between 
50 and 5000 MHz. They suggested that the lack of any radius- 
to-frequency mapping could be due to the emission originating 
from an 'inner cone', which is located closer to the 'core' emis- 
sion component, and so the effects of radius-to-frequency map- 
ping are less pronounced. 

In our model (see Figure|3) the two compo nents move close r 
together at low frequencies (in agreement wifh lLvne et al.l 19711) . 
We did not, however, see the co mponent separation increasing 
below 150 MHz (as reported bv ISieber et~aDll975l) . This be- 
haviour is similar to that of PSR B0809+74 at high frequencies, 
although in this case the components are only expected to pass 
through each other at ~ 1 MHz. Another curious feature of our 
model is that whilst component 2 gets broader at low frequencies 
(as expected), component 1 appears to get narrower. Radius-to- 
frequency mapping suggests that both components should get 
broader at low frequencies, as the emission region gets wider. 
In an attempt to conform with this idea, we tried to fit a three- 
component model. This fit the pulse profile reasonably well, 
and also agreed better with the radius-to-frequency model; the 
components didn't move closer together, and stayed at a con- 
stant width. However, a suitable fiducial point could not be de- 
termined, and there were large systematic errors in our timing 
residuals so the model was rejected. 

8. Discussion 

8. 1 . Profile Evolution 

We have shown that pulse profile evolution can introduce large 
errors into pulsar timing data, in agreement with the work of 
Ahuja et al] d2007h . These errors can be as large as a few mil- 
liseconds in some cases, depending on the period of the pulsar, 
and how asymmetric the pulse profile is. A frequency-dependent 
model of the pulse profile can be used to reduce these errors. 
Using this method, it was possible to define an analytical fidu- 
cial point in the pulse profile of each of the pulsars in our sample. 
This fiducial point is valid to within a few milliseconds (corre- 
sponding to ~ 1 degree in pulse phase), although the model does 
not remove the timing errors completely. Small timing errors re- 
main because the model is not an exact fit to the observed profile, 
and subtle differences between the shape of the modelled tem- 
plates and the data lead to systematics in the cross-correlation. 
Further work needs to be done to explore how to better remove 
these timing errors, and how they vary with time. 

We have found that radio emission from all of the pulsars 
in our sample originates from a narrow range of heights in the 
magnetosphere. The narrow ranges found do not fit well with 
models of radius-to-frequency mapping in a dipolar magnetic 
field. In addition, all of the pulsars in our sample show at least 
one other trait which cannot be explained by radius-to-frequency 
mapping. 

The asymmetric cones which we observe in PSR B0329+54 
were also observed by GG, who attributed their asymmetry to 
aberration and retardation. However, we have not detected any 
aberration or retardation effects in our timing residuals and we 
also find that the emission from the inner cone seems to be 
concentrated to within 183 km of the neutron star surface. The 
(~ 5 degree) skew in the cone corresponds to a time difference 
of ~ 10 ms, which is much greater than the aberration and re- 



Table 7. Parameters of the models used in the dynamic tem- 
plates. The functions given are appropriate for frequencies mea- 
sured in MHz, and pulse phase measured in degrees. Component 
brightnesses are defined relative to a fiducial component. 
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tardation effects possible from within this height range, which 
are ~ 1 ms (from Equation [T9|>. The fact that the outer cone is 
skewed in the opposite direction to the inner cone, also suggests 
that this cannot be explained by the standard model of a pulsar. 

PSR B0809+74 has a component which starts out as a pre- 
cursor at high frequencies and then drifts through the centre of 
the pulse profile, swapping sides with the central component and 
appearing as a postcursor at low frequencies. The frequency de- 
pendence of the position of the drifting component suggests that 
either refraction or some relationship between frequency and 
height (a change in height could explain a component being de- 
layed) significantly influences one component, but is not seen in 
the other. 

PSR B 1 13 3+16 shows emis s ion fro m a very narrow range of 
heights and, as IMitra & Rankinl d2002l) showed, component sep- 
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aration which increases too rapidly to be produced by dipolar 
field lines. One explanation for this could be that there are other 
mechanisms at work, which act together with the traditional pic- 
ture of a pulsar, complicating pulse profile evolution. 

PSR B 19 19+21 has a profile whose width decreases at lower 
frequencies. This is the exact opposite of what is predicted by 
radius-to-frequency mapping, and so is very difficult to explain 
using the standard picture of a pulsar. Again, there is a clear 
relationship between pulse shape and frequency, but it does not 
seem to be explainable by radius-to-frequency mapping. 

The fact that none of the pulsars in our sample behave as 
predicted by radius-to-frequency mapping suggests that a more 
complicated model of the pulsar magnetosphere is needed to 
describe pulse profile evolution. Although radius-to-frequency 
mapping has been successful in explaining some of the features 
seen in pulse profiles, it is clear that it cannot be used to fully 
describe any of the pulsars in our sample. There are, however, 
alternative theories which could potentially provide good fits to 
observational data. 



In the family of models proposed by Barnard & Aronsl 



(198 6 p, and de velope d furt her by Petroval 
Welt evrede etall (120031) and iBeskin & Philippovl (1201 ll) . 
the frequency-dependent profile evolution seen in pulsars is 
explained by propagation effects in the pulsar magnetosphere. 
In these models, the radiation originates from a small region in 
the magnetosphere, and refraction, dispersion and the different 
propagation modes (i.e. extraordinary and ordinary) in the 
magnetosphere are responsible for the frequency evolution of 
the different components which are observed in pulse profiles. 

iKaras tergiou & Johnston! d2007l) also provide an interesting 
empirical model of the pulsar magnetosphere, which could be 
used to explain all of the features which we have observed. They 
postulate that all radio emission originates from a patchy cone 
bounded by the last open field lines and that emission can come 
from any height, independently of frequency. Complex pulse 
profiles can then be explained by invoking emission from a range 
of heights, rather than assuming that the pulse profile probes the 
longitudinal shape of the beam at one single height. 

What causes pulse profile evolution is still an important 
question, and will be vital in understanding the pulsar emission 
mechanism, and for studies of pulsar geometries in the future. 
At this stage, it is still difficult to discriminate between the many 
models that exist, but next generation telescopes, like LOFAR, 
will be excellent tools for studying this effect. 



8.2. Magnetospheric Effects 

The argument that a more sophisticated model is needed to de- 
scribe radio emission from pulsars is also supported by consid- 
erations of aberration and retardation effects on our data. From 
these arguments, we have shown that radio emission from all 
of the pulsars in our sample is confined to a very small re- 
gion in the pulsar magn etos phere, which su pports the ideas of 
iBarnar d & Arons (1986) and lPetroval dlOOO). 



,,-2 



We have also shown that, as there is no departure from a 
dispersion law in our data, there is no evidence for super- 
dispersive delays or refraction from within the pulsar magne- 
tosphere in any of the pulsars in our sample. However, whilst 
we don't see a frequency-dependent delay in the timing residu- 
als, refraction may be needed to explain the broad component of 
PSR B0809+74, which is seen to drift through the pulse profile. 



8.3. ISM Effects 

In our data, we see no evidence for any deviation from the cold 
plasma dispersion law, suggesting that second-order ISM effects 
in these pulsars introduce additional time delays < 50 ns at nor- 
mal pulsar timing frequencies (1400 MHz). The parameters of 
the pulsars in our sample are typical of those found in the PTAs, 
the only difference being their longer pulse periods. This sug- 
gests that the ISM may not cause as much of a hindrance to 
pulsa r timing projects as first feared dHemberger & S tinebring 
I2008h . 

The fact that no unexpected delay was detected in any of 
our observations also means that (at least along these lines of 
sight) the ISM appears to be relatively smooth, with no large, 
dense structures. These findings agree with the idea that scatter- 
ing is dominated by one o r two small, but relativ ely high den- 
sity re gions as discussed by Stinebring (2006) and Brisken et alj 
d2010h . We have determined an upper limit on ^DM, which 
can be used as an indicator as to how much the DM is likely 
to change in the future. Comparing these predictions with real- 
ity will be a useful cross-check of how well this relation works. 
We have also been able to place upper limits on the scattering 
measure and the emission m easure. Although the se limits are 
weak compared to NE2001 dCordes & Lazioll2002l) . pulsar tim- 
ing could provide an independent method of measuring both of 
these quantities in the future. 



8.4. Future Observations 

The constraints set in this paper will be improved significantly 
by taking similar observations when LOFAR is completed, us- 
ing more stations at a lower observing frequency. The sensitivity 
of LOFAR has already improved by a factor of five since the 
observations for this paper were taken, and it is expected that it 
will increase significantly again soon, when the core stations can 
all be combined coherently. Increased sensitivity, particularly in 
the LB As, will reduce the error bars seen at low frequencies in 
our timing residuals, which dominate the uncertainty in our mea- 
surement. 

We could also increase our precision by observing at lower 
frequencies. LOFAR will soon be able to routinely observe with 
high sensitivity at frequencies as low as 15 MHz (see Figure flTTi. 
where the second-order ISM delays are expected to be at least an 
order of magnitude larger than they are at 40 MHz. 

By repeating this experiment in the future on the same set of 
pulsars, we could test whether pulse profile evolution is stable 
with time, and also track variations in the DM with great accu- 
racy. Both of these parameters are not completely understood, 
and are vital for high-precision pulsar timing. It would also be 
interesting to perform this experiment on a millisecond pulsar. 
A faster rotation rate and a narrower pulse (in absolute terms) 
means that TOAs can be determined more accurately, which 
would improve our constraints by at least an order of magnitude. 
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Fig. 3. PSR B 19 19+21 is the simplest of our models, a double peaked pulsar. The top left panel shows the model used for the pulsar, 
two Gaussians with the fiducial point set as the midpoint between the peaks of each component. The Gaussians are fitted to the data 
using a least squares fitting algorithm, allowing their width, height and separation to vary. This fitting is shown in the right hand 
panel. The data is plotted in light grey, the two fitted components are plotted in dark grey, and the sum of both fitted components 
is plotted in black. The Gaussian parameters for each of the observations are recorded and plotted as a function of frequency in the 
bottom left panel. These parameters are then fitted with power laws to get a model of the pulse profile as a function of frequency. 
This global model is subsequently used to produce templates for cross-correlation. The 'ratio of peaks' plot has discontinuities 
because different power laws were fitted in each observing band. PSR B 19 19+21 was not detected in the Effelsberg observations as 
the source was too weak. 
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Fig. 4. A plot of frequency against apparent DM in each of the observing bands. The circluar points show the DM obtained using 
a static template, and the triangular points show the DM from a frequency-dependent template based on our models of the pulsars' 
pulse shape evolution. In each case (apart from that of PSR Bl 133+16) the model-based template provides a consistent DM across 
all frequency bands, while use of the static template often results in significantly different DM values. 
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Fig. 5. Simulations of structure in our residuals of PSR B 19 19+21. The top panel shows an example of ourTOAs, with a simulated 
v~ 4 ISM-like delay, which would be 0.84 ms at 48 MHz, added to them. The fit to the data is shown by the black line, and the null 
hypothesis (no ISM delay) is plotted in grey. Similarly, the right hand panel shows an example of a simulated aberration/retardation- 
like delay with a v° 6 dependence which is 0.28 ms at 180 MHz. Because the errors on the TOAs are much smaller at high frequencies, 
we are much more sensitive to delays at high frequencies, despite the fitted jump. Note that the sensitivity of both the LOFAR LBA 
and HBA observations is now vastly improved and so we should be able to better constrain (or even detect) some of these effects in 
the near future. 
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